#ifndef HLLC_SCHEME_NULLARY_HPP
#define HLLC_SCHEME_NULLARY_HPP

#include "bi_hllc.hpp"

namespace hms::swe::hllc
{

template<typename T>
using NestedOrScalar = std::conditional_t<
	isDenseType<remove_qualifiers<T>>(),
	typename remove_qualifiers<T>::Nested,
	scalar
>;

template<typename T>
constexpr bool scalarOr1xN(){
	using U = remove_qualifiers<T>;
	if constexpr ( isDenseType<U>() ){
		return U::RowsAtCompileTime == 1;
	} else {
		return std::is_same_v<scalar, U>;
	}
}
 
template<
	typename N,
	typename Dl,  typename Dr,
	typename Vnl, typename Vnr,
	typename Vpl, typename Vpr,
	typename Sl, typename Sr>
class FluxFunctor
{
public:
	using Scalar = scalar;

	static constexpr Index
		RowsAtCompileTime { 3 },
		MaxRowsAtCompileTime { 3 },
		ColsAtCompileTime { N::ColsAtCompileTime },
		MaxColsAtCompileTime { N::MaxColsAtCompileTime };

	using DenseType = Array<
		Scalar, RowsAtCompileTime, ColsAtCompileTime,
		Eigen::ColMajor, MaxRowsAtCompileTime, MaxColsAtCompileTime
	>;

private:
	static_assert( isDenseType<N>() );

	const NestedOrScalar<N> m_normal;
	const NestedOrScalar<Dl> m_dL;
	const NestedOrScalar<Dr> m_dR;
	const NestedOrScalar<Vnl> m_vnL;
	const NestedOrScalar<Vnr> m_vnR;
	const NestedOrScalar<Vpl> m_vpL;
	const NestedOrScalar<Vpr> m_vpR;
	const NestedOrScalar<Sl> m_sL;
	const NestedOrScalar<Sr> m_sR;

	static_assert(
		N::RowsAtCompileTime == 2
		&& scalarOr1xN<Dl >() && scalarOr1xN<Dr>()
		&& scalarOr1xN<Vnl>() && scalarOr1xN<Vnr>()
		&& scalarOr1xN<Vpl>() && scalarOr1xN<Vpr>()
		&& scalarOr1xN<Sl >() && scalarOr1xN<Sr>()
	);

	using S = WaveRegion;

public:
	FluxFunctor(
		const N& normal,
		const Dl& dL, const Dr& dR,
		const Vnl& vnL, const Vnr& vnR,
		const Vpl& vpL, const Vpr& vpR,
		const Sl& sL, const Sr& sR
	) :
		m_normal { normal },
		m_dL { dL },
		m_dR { dR },
		m_vnL { vnL },
		m_vnR { vnR },
		m_vpL { vpL },
		m_vpR { vpR },
		m_sL { sL },
		m_sR { sR }
	{
		assert1xN_or_arithmetic( m_dL, m_normal.cols() );
		assert1xN_or_arithmetic( m_dR, m_normal.cols() );
		assert1xN_or_arithmetic( m_vnL, m_normal.cols() );
		assert1xN_or_arithmetic( m_vnR, m_normal.cols() );
		assert1xN_or_arithmetic( m_vpL, m_normal.cols() );
		assert1xN_or_arithmetic( m_vpR, m_normal.cols() );
		assert1xN_or_arithmetic( m_sL, m_normal.cols() );
		assert1xN_or_arithmetic( m_sR, m_normal.cols() );
	}

	template<WaveRegion side>
	Scalar massFlux(Index col) const {
		return (side == S::left ?
			hllc::massFlux(m_dL[col], m_vnL[col] ) :
			hllc::massFlux(m_dR[col], m_vnR[col] ) );
	}

	template<WaveRegion side>
	Scalar momentumFlux(Index col) const {
		return (side == S::left ?
			momentumFlux_normal(m_dL[col], m_vnL[col] ) :
			momentumFlux_normal(m_dR[col], m_vnR[col] ) );
	}

	template<WaveRegion side>
	Scalar momentum_norm(const Scalar& f1, const Scalar& f2, Index col) const {
		return m_normal(0, col) * f2 - m_normal(1, col) * f1 *
		( side == S::left ? m_vpL[col] : m_vpR[col] );
	}

	template<WaveRegion side>
	Scalar momentum_para(const Scalar& f1, const Scalar& f2, Index col) const {
		return m_normal(1, col) * f2 + m_normal(0, col) * f1 *
		( side == S::left ? m_vpL[col] : m_vpR[col] );
	}

	template<WaveRegion side>
	Scalar hllMassFlux(Index col) const {
		if constexpr ( side == S::left ){
			return hllFlux(
				m_sL[col], m_sR[col],
				massFlux<S::left>(col),
				massFlux<S::right>(col),
				m_dL[col], m_dR[col]
			);
		} else {
			return hllFlux(
				m_sR[col], m_sL[col],
				massFlux<S::right>(col),
				massFlux<S::left>(col),
				m_dR[col], m_dL[col]
			);
		}
	}

	template<WaveRegion side>
	Scalar hllMomentumFlux(Index col) const {
		if constexpr ( side == S::left ){
			return hllFlux(
				m_sL[col], m_sR[col],
				momentumFlux_normal(m_dL[col], m_vnL[col]),
				momentumFlux_normal(m_dR[col], m_vnR[col]),
				massFlux<S::left>(col),
				massFlux<S::right>(col)
			);
		} else {
			return hllFlux(
				m_sR[col], m_sL[col],
				momentumFlux_normal(m_dR[col], m_vnR[col]),
				momentumFlux_normal(m_dL[col], m_vnL[col]),
				massFlux<S::right>(col),
				massFlux<S::left>(col)
			);
		}
	}

	template<WaveRegion side>
	Scalar starFlux(Index row, Index col) const {
		switch(row){
			case 0:
				return hllMassFlux<side>(col);
			case 1:
				return momentum_norm<side>(
					hllMassFlux<side>(col),
					hllMomentumFlux<side>(col),
					col
				);
			case 2:
				return momentum_para<side>(
					hllMassFlux<side>(col),
					hllMomentumFlux<side>(col),
					col
				);
			default: UNREACHABLE;
		}
	}

	template<WaveRegion side>
	Scalar outerFlux(Index row, Index col) const {
		switch(row){
			case 0:
				return massFlux<side>(col);
			case 1:
				return momentum_norm<side>(
					massFlux<side>(col),
					momentumFlux<side>(col),
					col
				);
			case 2:
				return momentum_para<side>(
					massFlux<side>(col),
					momentumFlux<side>(col),
					col
				);
			default: UNREACHABLE;
		}
	}

	bool useFluxLeft(Index col) const {
		scalar
			helpL { m_dL[col] * ( m_vnL[col] - m_sL[col] ) },
			helpR { m_dR[col] * ( m_vnR[col] - m_sR[col] ) };
		return 
			sign(m_sL[col] * helpR - m_sR[col] * helpL) ==
			sign(helpR - helpL);
	}

	Scalar operator() (Index row, Index col) const {
		if ( m_sL[col] < 0 && m_sR[col] > 0){
			if ( useFluxLeft(col) ){
				return starFlux<S::left >(row, col);
			} else {
				return starFlux<S::right>(row, col);
			}
		} else if ( m_sL[col] >= 0 ){
			return outerFlux<S::left >(row, col);
		} else {
			return outerFlux<S::right>(row, col);
		}
	}
};

template<typename ... Ts>
using FluxReturnType = Eigen::CwiseNullaryOp<
	FluxFunctor<Ts ...>,
	typename FluxFunctor<Ts ...>::DenseType
>;

template<typename N, typename ... Ts>
auto flux( const Eigen::DenseBase<N>& normal, Ts&& ... args )
	-> FluxReturnType<N, Ts ...>
{
	using DenseType = typename FluxFunctor<N, Ts ...>::DenseType;
	return DenseType::NullaryExpr(
		3, normal.cols(),
		FluxFunctor<N, Ts ...>( normal.derived(), std::forward<Ts>(args) ... )
	);
}

} // namespace hms::hllc

#endif
